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Abstract. In this article, we develop goal-oriented error indicators to drive adaptive 
refinement algorithms for the Poisson-Boltzmann equation. Empirical results for the sol- 
vation free energy linear functional demonstrate that goal-oriented indicators are not suf- 
ficient on their own to lead to a superior refinement algorithm. To remedy this, we pro- 
pose a problem-specific marking strategy using the solvation free energy computed from 
the solution of the linear regularized Poisson-Boltzmann equation. The convergence of 
the solvation free energy using this marking strategy, combined with goal-oriented re- 
finement, compares favorably to adaptive methods using an energy-based error indicator 
Due to the use of adaptive mesh refinement, it is critical to use multilevel preconditioning 
in order to maintain optimal computational complexity. We use variants of the classi- 
cal multigrid method, which can be viewed as generalizations of the hierarchical basis 
multigrid and Bramble-Pasciak-Xu (BPX) preconditioners. 
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1. Introduction 

The Poisson-Boltzmann equation (PBE) is a widely used model for electrostatic inter- 
actions of charged bodies in dielectric media, such as molecules, ions, and colloids, and 
thus is of importance in many areas of science and engineering, including biochemistry, 
biophysics, and medicine. (See the classical texts [65, 74] for a derivation of the PBE.) 
The importance of the PBE model is reflected by the popularity of software packages 
such as APBS [9], CHARMM [34], DelPhi [70], and UHBD [64], within the molecular 
modeling communities. It provides a high fidelity mean-field description of electrostatic 
interactions and ionic distributions of solvated biomolecular systems in equilibrium. The 
partial differential equation itself is challenging to solve numerically due to singulari- 
ties of different orders at the positions of permanent point charges and the presence of a 
dielectric interface. 

In this article, we develop an adaptive multilevel finite element method for the PBE 
using goal-oriented a posteriori error indicators. This adaptive algorithm, which is a 
variant of that studied for the PBE in [55, 8, 9], deviates substantially from previous work 
in that the error indicator is based on a user defined quantity of interest or goal. This is 
in contrast to traditional residual-based adaptive refinement algorithms (like those devel- 
oped for the PBE in [39]) that drive -refinement to minimize the global error measured in 
an energy-norm. The goal-oriented refinement methodology has been successfully em- 
ployed in a wide range of application areas, including fluids, elasticity, and fluid structure 
interaction [11]. Despite these successes, we show that this methodology applied directly 
to the PBE does not necessarily lead to a successful adaptive algorithm. To remedy this 
issue we propose a novel marking strategy which recovers the performance commonly 
seen in other applications. This is the first time that this particular goal-oriented re- 
finement strategy has been applied to the PBE specifically, and molecular biophysics in 
general. 

At the core of any adaptive finite element approach are the iterative methods used to 
solve the discretized equation. However, due to the ill-conditioning of the linear systems 
arising from the discretization of the PBE, the convergence rate of traditional iterative 
solvers is significantly deteriorated. To remedy this, we combine modem Bramble- 
Pasciak-Xu (BPX)-type multilevel preconditioners with the goal-oriented adaptive al- 
gorithm mentioned above. When applied to the PBE, our results demonstrate that the 
overall algorithm is accurate, highly efficient and scalable with respect to the number of 
levels in the adaptive hierarchy. 

An outline of the article is as follows. In section 2, we give a brief overview of the 
Poisson-Boltzmann equation, and describe the most useful formulations for modeling 
and numerical simulation, such as the regularized formulations described in [39, 56, 37]. 
We also discuss the solvation free energy functional corresponding to a given reaction po- 
tential, which will form the basis of our goal-oriented error indicators developed later in 
the article. We describe adaptive finite element methods in section 3, including weak for- 
mulation of the regularized PBE, discretization by finite element methods, and adaptive 
algorithms driven by a posteriori error indicators. In section 4, we describe a particular 
class of error indicators known as goal-oriented indicators, and describe several indi- 
cators designed for the PBE. In section 5, we discuss a local multigrid algorithm used 
to precondition an iterative Krylov method for solving the linear systems arising from 
adaptive mesh refinement. With some care, these methods enable an algorithm whose 
complexity is close to optimal. The results from a sequence of numerical experiments 
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Figure 1 . Debye-HUckel model of a charged biological structure im- 
mersed in a solvent containing mobile ions 

using the Finite Element ToolKit (FETK) are presented in section 6. These results high- 
light the efficacy of the goal-oriented error indicator for the Poisson-Boltzmann problem, 
as well as the utility of the linear solver strategy combined with the adaptive algorithm, 
driven by the goal-oriented error indicator. We draw some conclusions in section 7. 



2. The Poisson-Boltzmann Equation 

The Poisson-Boltzmann equation (PBE) is a second-order nonlinear partial differen- 
tial equation whose solution gives the electrostatic potential, for a solute molecule 
immersed in an implicitly defined solvent. Using a mean-field approximation, the sol- 
vent is treated as a bulk medium where ions are distributed according to the Boltzmann 
distribution. Figure 1 is a schematic representation of the domain of the PBE, denoted by 
Vt. The innermost region, Vtm, contains the explicitly represented solute molecule. The 
outer region, Vtg, is the bulk solvent and contains the implicit solvent ions. Between Vtg 
and Vtm is the ion exclusion layer, which separates the solute from the solvent ions, and 
has a width dependent on the size of the solvent ions. For simplicity, we will assume the 
solvent ions are small and the ion exclusion layer can be neglected. Hence, the interface 
between the solute and solvent is a surface, denoted by F = 0^ fl Vtg- The shape of the 
surface is governed by the short-range repulsive van der Waals interactions, which pre- 
vent the solvent from penetrating the solute. The precise definition of the surface varies 
depending on the model [10]. 

The PBE for a 1:1 electrolyte (e.g., sodium chloride) is 



—V ■ e{x)Vu{x) + K^{x) sinh {u{x)) 



i=l 



e{x) 



u{oo) = 
du{x 
dn 



0, xeT, 



(2.1) 



where u{x) = ec(f){x) / is the dimensionless potential, Cc is the charge of an electron, 
ks is Boltzmann's constant, and T is the temperature. Here, [[ ■ ]] denotes the jump 
across the interface 



fix) 



lim /(x + Cn) - fix - Cn) 



(2.2) 
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and n is the outward pointing normal of dVLm- The dielectric function e(x) jumps one 
or two orders of magnitude at the interface V. For example, commonly used values are 
e{VLm) = Cm = 2 and e(^7s) = Cs = 80. The modified Debye-Hiickel parameter, R, has 
a similar discontinuity, with R{^lm) = and /^(^7s) = Rg > 0. The fixed ions within the 
solute are represented by a sum of Dirac delta distributions, with fixed charge centers, 
Xi, and charges CcQi. This charge distribution induces singularities in the electrostatic 
potential and has, until recently, proved to be difficult to treat numerically. 

To address this issue the PBE is reformulated so that the singularities are explicitly 
removed [48, 86]. Following [39, 56, 37], this is accomplished by writing the potential 
as a sum of a singular term Uc and a nonsingular remainder Ur- The singular term is the 
Coulomb potential 

1=1 ' 

which satisfies the Poisson equation 

Ane ^ 

—V ■ emVuJx) = - — > QiSix — Xi) for x e fl, 

^ ^ keT ^ ^ ' (2.4) 

i=l 

Mc(oo) = 0. 

There are numerous fast algorithms (with linear or near linear complexity) for evaluating 
Uc on a set of quadrature points, such as fast multipole [49], multilevel summation [72, 
51] and particle mesh Ewald [46]. 

Substituting u = Ur ^ Uc into the PBE (Eq. 2.1) gives a modified form of the PBE, 
which was termed in [39] as the Regularized PBE (RPBE): 

—V ■ e(a;)VMr(a;)+K^(x) sinh {ur{x) + Uc(x)) 

= V ■ (e(x) - em)VMc(x), xefi^Ufi^, 
Ur{po) = 0, (2.5) 



dUr{x) 



xer. 

on 



Note that, because both R{x) and (e(x) — are zero for x G Vim and the centers of 
the atoms in the solute are well separated from F, the singular function, Uc, is never 
evaluated near the singularities. This formulation was used in [39] to develop continuous 
a priori L°° estimates of solutions, and subsequently to show existence and uniqueness of 
solutions of the PBE. Furthermore, the authors established discrete a priori L°° estimates 
for Galerkin solutions, making possible quasi-optimal a priori error estimates, as well as 
a provably convergent adaptive finite element method for the RPBE. 

Recently, an alternative 3-term splitting of the PBE has been proposed in [56] which 
addresses the inherit subtractive cancellation in the reconstruction of the electrostatic po- 
tential used by the RPBE. In this article, the authors establish mathematical results for the 
alternative splitting, including continuous and discrete a priori estimates, existence 
and uniqueness of solutions, quasi-optimal a priori error estimates, and a convergent 
adaptive finite element method (AFEM). (Whereas in [39] only AFEM convergence was 
shown, it was shown in [56] that AFEM is a contraction for the RPBE, using a new 
AFEM convergence framework for nonlinear problems developed in [57].) The 3-term 
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splitting decomposes the electrostatic potential into 

where Uc is the Coulomb potential, but here it is restricted to the subdomain Vim- The 
harmonic term, Uh, is defined as the solution to 

-V^Uh{x) = inQrn (2.7) 

Uh{x) = —Uc{x) onF. (2.8) 

Applying the definitions of u^, Uh and substituting into the PBE (Eq. 2.1), one obtains an 
equation for u^, 

-V ■ {e{x)Vu^{x)) + K^{x) sinh(M3(x)) = in (2.9) 



■ .duz{x) 



d{uc{x) + Uh{x)) 

em 7. onT, (2.10) 

on 



M3(CX))=0. (2.11) 

In contrast to the RPBE, this formulation avoids the subtractive cancellation since = 
u in Vtg, and the Coulomb term is not used to reconstruct the potential in the solvent 
subdomain. 

For systems which are not highly charged, the variation in the potential is relatively 
small, and the hyperbolic sine term is well approximated by its linearization. This 
approximation, which replaces sinh(M) with u, results in what is known as the lin- 
ear Poisson-Boltzmann equation (LPBE), or the linear regularized Poisson-Boltzmann 
equation (LRPBE) in case of the RPBE. Although this approximation reduces the ionic 
response of the solvent [10, 47], it can significantly reduce the complexity of many nu- 
merical algorithms (e.g., boundary element and boundary integral methods) [63, 25, 62]. 

One important use for the solution to the PBE is in the calculation of solvation free en- 
ergies. This quantity measures the thermodynamic work of moving the solute molecule 
from a vacuum to a solvent environment. The solvation free energy can be written as a 
sum of nonpolar and polar contributions. The nonpolar term depends on the solvent ac- 
cessible surface area, excluded volume, and nonpolar forces which are typically assumed 
to be independent of the electrostatic potential [60, 78]. The polar term, S, is a linear 
functional of the solution to the RPBE, Ur, (also known as the reaction potential) [60] 
and can be expressed as 



S{Ur) 2 



1 r ^ 

- I Ur{x)y^ecqi8{x - Xi) dx. (2.12) 



For the 3-term splitting, the reaction potential is the sum of and Uh- 



3. Adaptive Finite Element Methods 

The finite element approach provides a natural framework for dealing with the com- 
plex molecular surfaces which arise in the PBE. Although there are modified finite dif- 
ference methods which address this difficulty [83], finite element methods provide an 
attractive alternative when paired with an adaptive unstructured mesh designed to con- 
form to the shape of the solute molecule [62]. In this section, we present a general 
adaptive finite element method for the regularized PBE, including the weak formulation, 
discretization, solution using an inexact global Newton iteration, and adaptive refinement 
procedure. For more details on the finite element method, see [73, 32, 26, 45]. 
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3.1. Weak Forms. To give a well-defined weak formulation, the nonlinearity involving 
exponentials must be controlled; in [39, 56], a priori estimates are obtained for any 
solution to the RPBE, giving almost everywhere pointwise bounds of the form: a ^ 
Ur ^ (3. This leads to working with a well-defined solution space that consists of a 
non-empty, topologically closed, convex subset of H^{^1): 

Me := {v E H^{n) : a ^ v ^ /3 a.e. mn,v = u - UcOn dn}. (3.1) 

It is shown in [39, 56] that there exists a unique solution to either regularized form of the 
RPBE in Mg C if^(fi). The weak formulation is: Find Ur E such that 

a{ur^v) + h{ur + UciV) = L{v) "iv E HliVt) (3.2) 

where 

a{u,v) = {eWu.Wv) (3.3) 

f ) = (k^ sinh(M), t>) (3.4) 

L{v) = {-{e~e^)Vu,,Vv). (3.5) 

The linear functional L{-)h defined by integrating the right hand side of Eq. 2.5 by parts 
and applying the jump condition to eliminate the interface terms. 

The weak form for the 3-term split regularized PEE requires solving two problems: 
first for the harmonic term on VL^ and second for the split potential on the whole domain 
VL. Define the solution space to the harmonic problem as := {v E H^{VLrn) '■ v{x) = 
—Uc{x) Vx E d^lm}- Then the weak form of Eqns. 2.7-2.11 is: Find (uh, E x Me 
such that 

a{us,v) + b{u3,v) + am{uh,w) = {g{uh),v) \/{w,v) E H^in^) y< H^{^) (3.6) 
where am(-, ■) is the restriction of the bilinear form to the subdomain and 

{g{uh),v) = I e^ ^^^''^^''K dx. (3.7) 
Jdn,„ on 

3.2. Solving. Due to the hyperbolic sine, the RPBE has a strong nonlinearity. The dis- 
cretized nonlinear problems defined in Eqns. 3.2 and 3.6 can be solved using an inexact- 
Newton method [44]. For brevity, we give details for Eq. 3.2. Define the weak residual 
functional to be 

{Riu'i^^v) = L{v) - (a(M^ v) + + M^, v)) . (3.8) 
Here is the discrete solution satisfying the system of nonlinear equations 

{R{u^^),v) = \/veV'' (3.9) 

where is the space of piecewise linear functions defined by the tetrahedral mesh. 
Linearizing Eq. 3.9 around u'^ results in 

Jw'' := {DR{u'^)w\v) 

= -a{w\v) -b'iu';: + Uc;w\v) VveV''. (3.10) 

In the linear RPBE, sinh(n) is replaced by u, and b'{u, v) = {k^u, v). Newton's method 
defines the nonlinear update vector as the solution to 

Js^ = -R{u^^). (3.11) 
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Given an initial guess ^ uq, the updated solution is defined as ui = uq + Sh- This 
process can be repeated until a desired level of convergence is achieved. An inexact- 
Newton method uses an iterative solve to find an approximate solution to Eq. 3.1 1, with 
a relatively large tolerance for the linear solve when far from the nonlinear solution. 
However, as the exact solution to Eq. 3.9 is approached, the linear solver tolerance is 
tightened so that quadratic convergence is achieved. 

The computational complexity of the Newton solver is dominated by the method used 
to solve the N linear algebraic equations [15, 50] within each iteration. Multilevel meth- 
ods provide an advantage in that they are provably optimal or nearly optimal methods 
for solving these systems [13, 50, 80]. The presence of geometrically complex discon- 
tinuities in the dielectric e and in the Debye-Hiickel parameter R in the PBE destroy 
classical multilevel method efficiency, and can even cause divergence. This is analyzed 
at length for the PBE in [53, 58, 59], and various techniques based on coefficient av- 
eraging and algebraic enforcement of variational (Galerkin) conditions are examined. 
Algebraic multilevel methods have been used successfully for many similar problems; 
cf. [16, 17, 30, 36, 35, 31, 71, 76, 77]. A fully unstructured algebraic multilevel ap- 
proach is taken in FETK, more details are provided in Section 5. 

Starting from an initial mesh %, the adaptive mesh refinement procedure builds a se- 
quence of conforming meshes To, 7i, . . . , 71 [39, 56, 43]. This procedure is divided into 
four steps: SOLVE, ESTIMATE, MARK, and REFINE. In the SOLVE step, a solution is 
computed on the current mesh. Using this result, the ESTIMATE step computes elemen- 
twise error indicators and an estimate of the global error. In a production environment, 
the procedure terminates if the global error estimate is below some prescribed tolerance. 
Otherwise, the MARK step selects elements for refinement. This step is crucial to the 
convergence of the method. Finally, REFINE subdivides the marked elements possibly 
subdividing additional unmarked elements in order to produce a conforming mesh. The 
refinement technique used in this article is longest edge bisection [69]. 

4. Error Indicators 

In this section we present a posteriori error indicators for use in adaptive refinement. 
These estimators are typically developed by considering the residual of the weak form. 
For example, given a finite element solution E the weak residual for the linear 
RPBE is (compare to Eq. 3.8) 

v) = L{v) - (/cV, v) - (a(u^ v) + {k\^, v)) (4.1) 

for a given v e V. For the remainder of this article, we restrict our attention to two 
classes of error indicators: energy-based and goal-oriented. The first class estimates the 
error in the energy norm, although this idea can be generalized to other norms, (e.g., the 
norm). The second class, called goal-oriented, focuses on estimating the error in a 
user specified quantity of interest or goal functional. In the following sections, we derive 
error indicators from both classes for the linear RPBE. For the goal-oriented indicator, 
the solvation free energy is used as the target functional. 

4. 1 . Energy Norm Indicators. A standard a posteriori error indicator is based on bound- 
ing the error in the energy norm. It is easily derived by breaking the weak residual into its 
elementwise components and integrating by parts over each element [1]. This technique 
was used in [40] to derive the following estimator for the linear RPBE 

7]l{u'') = h^WrKWhf^K) + \heK\\rdK\\l2(QK), (4-2) 
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where 

rxix) =(V ■ (e(x) - em)Vuc{x) - r'^{x)Uc{x)) 

- {-V ■ e{x)Vu'^{x) + k\x)4{x)) ^KeT 

and 



(4.3) 



raxix) = riK ■ 



(e(x) - emWuJx) + e(x)Vu:(x) 



e r. (4.4) 



2 



This indicator gives a bound on the error measured in the energy norm (|||t;||| := a{v, v) + 
{k^v, v)) 

\\\Ur~u'r\f <Y.'^Ki^^r)- (4-5) 

KeT 

Bounding the error in other norms is possible. For example, in [39, 56] a similar a 
posteriori error indicator for the RPBE was shown to bound the error measured in the 
norm. Other efforts have focused on formulating the RPBE as a first-order system 
least squares (FOSLS) problem, which has a natural error estimate [22, 38]. 

4.2. Goal- Oriented Indicators. Key to the development of goal-oriented error indica- 
tors is relating the weak residual to the error in the goal functional. For symmetric linear 
problems, a direct application of the Riesz representation theorem shows that there exists 
a dual function that when paired with the weak residual gives the error in the goal [21]. 
The challenge is to approximate this function and utilize that approximation to develop 
error indicators. However, for nonlinear problems, like Eq. 2.5, the definition of the dual 
function is not so clear. The first part of this section discusses a strategy for defining the 
dual function for both the nonlinear RPBE and the three term splitting. Using this def- 
inition of the dual, two goal-oriented error indicators are proposed for the linear RPBE 
utilizing the solvation free energy as the quantity of interest. 

Following [21], a dual function for the nonlinear RPBE can be defined by considering 
the constrained minimization problem 

= argmin ^(u*) subjectto a{ur,v) + b{ur + Uc,v) = L{v) e Hl{^), (4.6) 

where a(-, ■), b{-, •) and L(-) are specified in Eq. 3.5, and S is the goal functional. Note 
that because the RPBE constraint determines the solution uniquely, the minimization 
problem has the same solution. However, specifying the minimization provides an ad- 
ditional mathematical framework to define the dual function. To see this, consider the 
Lagrangian associated with the minimization problem 

Q{Ur, w) = S{Ur) + {L{w) — {a{Ur,w) + b{Ur + Uc, w))) , (4.7) 

where the Lagrange multiplier, w E Hq^VI), is also the dual function. Taking the first 
variation of 9 with respect to u gives the dual problem: 

Find w e H^in), such that {DR{ur)v, w) = -S{v), E H^in), (4.8) 

where {DR(-)-, ■) was defined in Eq. 3.10. As discussed above, if b{-, ■) is linear in the 
first argument then Eq. 4.8 simplifies to a{v,w) + {r'^v,w) = S{v) E Hq^VI). For 
the linear problem the error in a goal functional S{-), like the solvation free energy, is 
simply expressed in terms of the weak residual 

S{Ur—U^.) = a{Ur — u'^,w) + {R'^{Ur — U^.),W) 

= L{w)-a{u^;,w)-{R\u^; + Uc),w). (4.9) 
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Thus, if w is known, the error in S{-) is easily calculated. In the nonlinear case the error 
in the goal satisfies 

S{ur - = L{w)- w) - + Uc] ul,w) + E, (4. 10) 

where E is quadratic in the error in [21]. 

For the three term splitting of the PBE, again we setup a constrained minimization 
problem to define the dual. Using the notation from the weak form in Eq. 3.6, the corre- 
sponding Lagrangian is 

^3-tena{i^3, Uh]W3,Wh) = Siu^) + Sm{uh) + 

{g{uh),W3) - (a(M3, ws) + b{u3,W3) + am{uh,Wh)) (4.11) 

where u'3 G Hq(^1) and Wh E ifQ(fim) are dual functions. The functional Smi-) is a 
restriction of the original goal functional to the f2m domain. Taking the first variation of 
Qs-term with respect to and gives the dual problem 

a{v,ws) + b'{us;v,ws) = S{v) \/v e H^{^), 

amiv, Wh) = Sm{v) + {g'iUh, V),W3) WV e HQ^Qrn) 

where 

dv 



(4.12) 



{g'{uh]v),W3) = em^wsdx. (4.13) 
Jan^ on 

4.2.1. Goal-Oriented Error Indicators for the Linear RPBE. To make the application 
of the dual functions in error indicators more concrete, we present two goal-oriented 
indicators for the linearized RPBE. As an example we will focus on accurate computation 
of the solvation free energy (see Eq. 2.12). Unfortunately, S{-) is not bounded on Hq(^1) 
due to the inclusion of delta distributions. A common approach to circumventing this 
issue is to use a mollified version of the functional [1, 11]. In this case, the mollified 
solvation free energy is 

S{ur) ^ S'^i^Ur) = - / Ur{x)''^^ecqi9{\x — Xi\,a) dx, (4.14) 
^ i=i 

where 6* is a locally supported function defined such that 

limj9{\x\,a)f{x)dx = J 5{x)f{x) dx = f{0). (4.15) 
One possible choice for 6 is the step function 

where B^^ is the volume of a ball of radius a. 

A simple error indicator suggested by Eq. 4.9 is to first solve the dual problem using 
the same approximation space as the primal. This approximate dual could then be sub- 
stituted for w in Eq. 4.9 to compute the value of the indicator. However, if the same finite 
dimensional space is used for solving both the dual problem and the primal problem, 
then, because of Galerkin orthogonality, Eq. 4.9 will be zero (see [11]). A remedy is to 
instead solve the dual problem using a finer approximation space, f/^ C Hq(^1). One 
convenient choice is to maintain the same mesh and use higher order polynomials for 
U^. In the examples below, V'^ is the space of piecewise linear polynomials and is 



10 B. AKSOYLU, S. BOND, E. CYR, AND M. HOLST 

the space of piecewise quadratics. Let the finer resolution solution of the dual problem 
be denoted w'*'^. Substituting w ~ w^'^ into Eq. 4.9 yields 

- {aiu'::,w'''^-Pl:^^w''^') + iR\u'^ + u,),w''^'-Pl:^^w''^^)) (4.17) 

Where Pll^2 ^ convenient projection (e.g., nodal injection) of the fine space U'^ onto 
V'^. The choice of the projection operator will affect the quality of the indicator. Decom- 
posing the error into its elementwise contributions gives 



K 

- {R'iu'^ + u,),w^^^ - PI,w^'^)k < J2vk{u',w'^''), (4.18) 

K 

where the subscript K indicates the restrictions of the linear functional, bilinear func- 
tional or inner product to element K and 



K 



e(x) - e™)VM,(x) ■ Wiw'^'^x) - w\x)) 



(4-19) 

-e{x)Vu':^{x)-V{w'''\x)-w''{x)) 



r\x)u^,{x){w''^\x) - w\x)) 



dx. 



Here, the absolute value of the integrand in Eq. 4.18 has been taken over each element. In 
the numerical experiments in section 6, Eq. 4.19 is referred to as the the "goal-quadratic" 
error estimator used by the various adaptive refinement marking strategies. 

The error indicator in Eq. 4.19 requires solving a dual problem which is substantially 
larger than the primal problem. To alleviate this issue, we develop a second goal-oriented 
error estimator that finds an approximation to the dual in V^' (which is the same space 
as the primal problem). The error in the goal is estimated by solving many local el- 
ementwise boundary value problems. The technique proposed here is similar to the 
development of the equilibrated residual method for computing goal-oriented estima- 
tors [66, 67]. However, the less accurate but simpler element residual method (ERM), as 
discussed in [1], is used. Using the parallelogram law, the error in the linear functional 
can be rewritten [1, 66, 67] as 

S''{Ur - U^.) = a{Ur - U^,W - W^) + {R^{Ur - U^),W - W^) 

1 ,1 , (4.20) 

= I IIIK - U'^r) + - ^'')||r - 4 IIIK - 4) -{W- W^)f 

where again = a(f , v) + {R^v, v) is the square of the energy norm. Define element- 
wise error functions (pK = Ur — u'^\K and iI^k = w — w'^ | x to be computed by the element 
residual method. The error in the primal problem on element K is approximated by the 
solution to 



= {R{v!::),v)k+ ! rKv{x)ds yveH\K), (4.21) 

JdK 

where 

{R{u';:),v)k = Lk{v) - aKiu';:, v) - {r\u';: + Uc),v)k, (4.22) 
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and 

= {e{x)Vu'^{x) + {e{x) - em)Vu,{x)) ■ uk- (4.23) 
Similarly, the error in the dual problem on element K is approximated by the solution to 

f (4 24) 

+ / f^v{x)ds Wv e H\K), 

where 

f^ = e{x)Vw\x)-nK. (4.25) 
For a derivation of these equations see [40]. Explicitly stated, the "goal-linear" error 
estimator used in the numerical experiments in section 6 is given by 

r^M, w^) = ^ Uk + i^Kll - \ Uk - ^kII ■ (4.26) 

5. Multilevel Preconditioning 

As the mesh is refined, the conditioning of the linear system deteriorates, and precondi- 
tioning is necessary to accelerate the convergence of iterative solvers (e.g., the conjugate 
gradient method). The challenge in designing an efficient preconditioner is balancing the 
cost of applying the preconditioner with its effectiveness in improving the conditioning 
of the underlying system. 

In a given finite element mesh, at level j, we denote the set of nodes and its cardinahty 
by J\fj and Nj, respectively. We call the set of nodes introduced precisely at level j 
the fine nodes, and denote them A/"/. As the mesh is refined, AfJ is appended to J\fj-i, 
leading to the following hierarchy of nodes: 

AT, =a/',_iUa/;^ j = i,...,j, 

where Nj = N. 

In the local mesh refinement setting, the way the coarse and fine nodes are processed 
plays a central role in determining the overall efficiency of a preconditioner. If the com- 
putational cost per level can be maintained proportional to Nj — Nj^i, or slightly larger, 
then total cost will be order N, and the resulting preconditioner is said to have opti- 
mal computational complexity per iteration. If the resulting preconditioned system has 
a bounded condition number (independent of problem size), a solution can be obtained 
using an iterative method with a bounded number of iterations. Hence, the combination 
of optimal per iteration complexity with a bounded condition number leads to a solver 
with optimal overall complexity. 

In this article, we restrict the presentation of local multilevel preconditioning to a 
purely geometric (node based) perspective because the computational complexity is ex- 
actly governed by the number of nodes processed by the preconditioner at each level. The 
local multilevel preconditioners of interest can be classified into two groups: multiplica- 
tive local multigrid (MG) [18, 20, 28, 68] and additive local MG [28, 29]. The additive 
local MG preconditioner is often called the Bramble-Pasciak-Xu (BPX) preconditioner 
in the literature. In this article, we report on only the multiplicative variants. We use the 
term "classical" to refer to the application of a preconditioner in the uniform refinement 
setting. There is abundant literature on MG preconditioners [33, 75, 80, 81] and local 
MG (e.g., see [19, 79] for a review). 

Proofs demonstrating the optimality of classical MG and classical additive MG pre- 
conditioners rely on a geometric increase in the number of nodes per level. This is 



12 



B. AKSOYLU, S. BOND, E. CYR, AND M. HOLST 



because the cost per iteration of these classical preconditioners is proportional to Nj (not 
Nj — Nj_i ) per level, which results in suboptimal complexity if Cj Nj is not 0{N). 
This occurs frequently in the local refinement setting due to the slow increase in the num- 
ber of nodes between levels. 

The family of hierarchical basis (HB) preconditioners, developed by Bank, Dupont, 
and Yserentant [12, 14, 82], maintain a per-level cost proportional to Nj — Nj^i, by only 
processing (smoothing) the fine nodes at each level. Although the cost per iteration is 
optimal, HB preconditioners do not achieve a uniformly bounded condition number, and 
suffer from (9( J^) and 0(2"') iterations in two- and three-dimensions, respectively. To 
address this deficiency, we investigate local MG preconditioners which process a larger 
set of nodes, Xj, but still maintain a cost which is proportional to Nj — iV^-i at each 
level. Hence, we seek a set Xj such that 

AfJ c X, c Af„ 

with cardinality, Xj, which is proportional to Nj — Nj^i. At the same time, Xj, should 
be large enough that the resulting system has a bounded condition number, leading to 
a solver with optimal overall complexity. Aksoylu and Hoist [6] showed that this is 
possible even for three-dimensional local refinement routines. 

In the local mesh refinement setting, Aksoylu, Bond, and Hoist [2] studied the imple- 
mentation and algebraic aspects (e.g., matrix representations) of multilevel precondition- 
ers. Subsequent articles provide a comprehensive overview of local MG preconditioners 
with various emphases: for a theoretical treatment, see [3, 4, 5]; for optimality analysis 
in three-dimensional local refinement routines, see [6]; for surface mesh applications in 
computer graphics, see [7]. 

5.1. Local multigrid preconditioners. As mentioned in the previous section, the fun- 
damental difference between classical and local MG preconditioners is the smoothing 
operation. In classical MG, the smoother acts on all degrees of freedom on every level. 
In contrast, local MG only smooths a small subset, typically a neighborhood, of the fine 
degrees of freedom. Pseudo-code for a local MG V-cycle is provided in Algorithm 5.1. 



coarsest level solve 

smooth Si times on Xj 
restrict residual 
set e[j_i] to zero 
coarse-grid recursion 
add interpolated correction 
smooth S2 times on Xj 

The set of nodes used by the underlying method's smoothing operation defines the type 
of preconditioner. In this article, we introduce three different sets of nodes to be used in 
the local MG preconditioner. On a given refinement level, the marked region is the set of 
elements which have been marked for refinement. The refinement region is union of the 
newly introduced elements as a result of the refinement and closure procedure. Figure 2 
depicts a two-dimensional refinement region formed by a local refinement routine which 



Algorithm 5.1. Local multigrid V-Cycle: 

U[j] = Vcycle{u[j]J[j]J) 

0) If 3 = 1. solve = f[j] 

and return My]/ 

1) My] = Smooth{uu] , , /yj , Xj , si); 

3) ey_i] = 0; 

4) ey_i] = Vcycle{e[j^i] , ry_i] , j - 1); 

5) My] = uy] + /gLi]ey_i]; 

6) uy] = Smooth{u[j] , Ay] , fy] , Xj , S2); 

7) return uu] 
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consists of quadrasection closed by bisection (the so-called red- green refinement). The 
sets of nodes used to define the preconditioners are as follows: 




HB: o. 

BPX: o, □. 

BEK: o, □, +. 

ONERING: o, □, +, A. 



Figure 2. An example of red-green refinement (quadrasection- 
bisection) in two dimensions. On the left, the mesh before refinement, 
with the marked region shaded in gray. On the right, the mesh after re- 
finement (red edges) and closure (green edges), with the refinement region 
shaded in gray. The labels indicate which nodes are assigned to each of 
the preconditioner smoothing sets, Xj. 

• Xj-HB: The set of fine nodes on level j, i.e., JVJ [12, 14, 82]. 

• Xj-BFX: The set of nodes whose corresponding basis functions have support 
entirely contained in the refinement region [28]. 

• A'j-BEK: The set of nodes whose corresponding basis functions have non-empty 
intersection with the marked region. This set is named after the Bornemann- 
Erdmann-Komhuber type refinement [23] routine. It consists of fine nodes and 
their immediate neighboring coarse nodes in the marked region. For PI elements, 
this set can be inferred from the nonzero pattern of the prolongation operator. 

• A'j-ONERING: The set of nodes whose corresponding basis functions have non- 
empty intersection with the refinement region [2, 24, 27, 41]. This set consists 
of fine nodes and their immediate neighboring coarse nodes in the refinement 
region. For PI elements, this set can be inferred from the nonzero pattern of the 
coarse-fine subblock of the stiffness matrix. 

As an example, we have labeled the nodes in Figure 2 to show which nodes are in each 
of the sets described above. The set corresponding to the classical multigrid precondi- 
tioner contains all of the nodes on each level. 

We should note that the practical implementation of the various local MG precon- 
ditioners varies significantly depending on the particular preconditioner. Special care 
must be taken in order to achieve optimal computational as well as storage complexities. 
The implementation aspects of how to construct optimal complexity preconditioners are 
studied in more detail in [2] . 

6. Numerical Experiments 

6.1. Adaptive Refinement. In this section, the effectiveness of goal-oriented mesh re- 
finement is compared to refinement using the energy-based error indicator. Of interest is 
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Figure 3. Convergence of the solvation free energy for Fasciculin-1 
with both goal-oriented and energy-based indicators using the global 
marking strategy. 




— Goal-Oriented 
• — • Energy-Based 



computing the solvation free energy of the 921-atom Fasciculin-1 protein [61], using the 
solution of the linear RPBE. All tests were performed using FETK [54]. 

In order to solve the RPBE, a definition of the molecular surface and a mesh con- 
forming to that surface is needed. Various definitions of the molecular surface have been 
proposed in the literature, and the particular value obtained for the solvation free energy 
will depend strongly on the surface geometry [10]. However, the performance of the 
algorithms proposed here are insensitive to the choice of molecular surface, as long as it 
is sufficiently smooth, and the underlying mesh is conforming. Historically, generation 
of the mesh conforming to the surface was a great impediment to using finite elements 
for solving the PEE, and only recently, with the development of tools like PDB2PQR 
and GAMer, has molecular meshing become a routine task. The first step is to prepare 
the structure using PDB2PQR [42], which adds missing hydrogens, assigns charges, and 
specifies a radius for each atom in the protein. Next, the resulting PQR file is passed to 
GAMer [56, 52, 84, 85], which produces a tetrahedral mesh conforming to the shape of 
the protein. Finally, this mesh is used by FETK to solve the linear RPBE, and compute 
the solvation free energy. 

6.1.1. Global Marking Strategy. The adaptive refinement algorithm creates a sequence 
of grids To, 7i, . . . 71 . . . based on an error indicator. Critical to this algorithm is the third 
step MARK. The goal of this step is to select elements for refinement. There are several 
different choices for marking strategies [1, 11]. The strategy used here is to mark all 
elements in the Z*'^ refinement level that satisfy 

t]k > 1 max r]T G 71, (6.1) 

where 7 G (0, 1). This criteria yields no refinement for 7 = 1 and uniform refinement 
for 7 = 0. 

The convergence of the solvation free energy for Fasciculin-1 using the global mark- 
ing strategy with both energy-based (Eq. (4.2)) and goal-oriented indicators (Eq. (4.19)) 
can be seen in Fig. 3. The figure shows the relative error in the solvation free energy as 
a function of the number of unknowns in the primal problem. The error in the solvation 
free energy is estimated by computing a high resolution solution to the PBE using uni- 
form mesh refinement. The vertical line near 1.5 x 10^ unknowns marks the size of this 
reference solution. Notice that this figure differs from traditional finite element conver- 
gence plots which show error as a function of element radius. Since the resolution of 
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Figure 4. A cut-away of the 3D mesh surrounding Fasciculin-1. The 
colors indicate the distribution of marked elements using the global mark- 
ing strategy with either the energy-based indicator (left) or the goal- 
oriented indicator (right). Red and blue are marked elements in the sol- 
vent and solute subdomains, respectively. 



an adaptively refined mesh can greatly vary over the spatial domain, the element radius 
is not an appropriate measure of the resolution of the mesh. Furthermore, the order of 
accuracy of the approximation is based on the type of basis functions used, and thus is 
the same for both uniform and adaptive refinement. The benefit of adaptive refinement 
is that the mesh can be refined in regions that heavily contribute to the error, resulting in 
higher accuracy with fewer total degrees of freedom. 

In Fig. 3, the lower line shows the convergence of the solvation free energy for several 
levels of adaptive refinement using the energy-based indicator with the global marking 
strategy. This scheme makes steady progress to the correct solvation free energy. On 
the other hand, the upper line shows the results using goal-oriented indicators with the 
global marking strategy. This scheme performs very poorly. 

The reason for the poor performance can be explained by looking at Fig. 4. The images 
are cut-aways of the 3D mesh, colored to indicate the distribution of marked elements. 
The left image shows elements selected by the global marking strategy using the energy- 
based indicator. While the image on the right uses the goal-oriented indicator. The white 
elements are unmarked elements in the solvent subdomain and the gray elements are 
unmarked elements in the solute subdomain. Elements colored red are marked solvent 
elements, while blue elements are marked solute elements. Notice for the energy-based 
indicator only elements in the solvent subdomain are marked. This indicates that the 
reaction potential in the solute domain is relatively well approximated compared to the 
solution in the solvent subdomain. The distribution of the elements marked using a goal- 
oriented indicator is focused on a few locations in the inner subdomain around the solute 
atoms. This explains the poor convergence for the goal-oriented indicator in Fig. 3. Since 
the strategy does not indicate there is error in the solvent domain, no refinement takes 
place. 

To gain further insight into why the global goal-oriented strategy only marks elements 
in the solute domain, we refer the reader to Fig. 5. This image shows a cut-away with 
elements colored by their approximate signed contribution to the error in the goal. This 
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Figure 5. A cut-away of the 3D mesh surrounding Fasciculin-1. The 
colors indicate the positive (green) and negative (orange) estimated ele- 
mentwise contributions to the error in the solvation free energy (see 
eq. (4.17)). 



is constructed using a quadratic approximation of the dual and the signed elementwise 
contributions from Eq. (4.17). Note that because these are error contributions and not in- 
dicators, they take both positive and negative values. The element contributions range in 
value between —0.6928 and 1.2779. In the figure, positive element contributions greater 
than 0.005 are colored green and negative contributions less than —0.005 are colored 
orange. From the image it is clear that the contributions in the solute subdomain have 
relatively large magnitude, and they oscillate in sign. The result of this oscillation is that 
the majority of these contributions cancel when integrated over the entire solute domain. 
However, the error estimator in Eq. (4.18) uses the absolute value, which results in an 
overestimation of the error attributed to the solute domain. As a result, the solute do- 
main is over refined, unless steps are taken to modify the marking (or error estimation) 
strategy. 

6.1.2. Split Marking Strategy. To improve on the convergence of goal-oriented refine- 
ment, a second marking strategy is employed. This is a domain dependent marking 
strategy that attempts to spread the refinement over solvent and solute regions of the do- 
main. The strategy relies on splitting the mesh into two subsets T{ d Ti and 7^™ C 71, 
where Tf and 7^"* contain the elements in the solvent and solute domains respectively. 
Stated concisely, the marking strategy is 

Mark all li^ such that ' (6.2) 

K G /; Vk > 1 max tit 

where 7 G (0, 1). 

The split-marking strategy marks a significantly different group of elements, especially 
for the goal-oriented indicators. The color coding of Fig. 6 is the same as in Fig. 4. 
Again the energy-based refinement selects elements primarily in the solvent subdomain. 
However, because split-marking forces refinement in both subdomains, a few elements 
along the interface in the solute domain are also selected. In contrast, split-marking 
using a goal-oriented indicator marks a few elements in the solvent subdomain, while 
also marking elements surrounding the solute atoms. 




Figure 6. A cut-away of the 3D mesh surrounding Fasciculin-1. The 
colors indicate the distribution of marked elements using the split marking 
strategy with either the energy-based indicator (left) or the goal-oriented 
indicator (right). Red and blue are marked elements in the solvent and 
solute subdomains, respectively. 
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Figure 7. Convergence of the solvation free energy for Fasciculin-1 
with both goal-oriented and energy-based indicators using the split mark- 
ing strategy. The convergence measured against the size of the pri- 
mal/dual problem is on the left/right. 

Figure 7 shows the relative error of the solvation free energy in Fasciculin-1 as a 
function of problem size for several indicators using the split-marking strategy. In the 
figure, there are two plots. The plot on the left shows the relative error in the solvation 
free energy as a function of the number of unknowns in the primal problem. The second 
plot, with the exception of the energy-based refinement strategy, shows the relative error 
as a function of the size of the dual problem. For energy -based refinement, since no dual 
is needed, the horizontal axis is the size of the primal problem. 

In the figure, the split marking strategy using the energy-based indicator converges 
steadily. This is similar to the convergence of energy-based refinement using global 
marking. On the other hand, there is a dramatic improvement in the convergence of 
both the linear and quadratic goal-oriented refinement techniques, with the rate even 
increasing slightly as the size of the problem increases. Compare this with the poor 
results for goal-oriented refinement with the global marking strategy from Fig. 3. The 
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Figure 8 . Number of primal unknowns as a function of the number of 
levels of refinement. The reference hne grows proportional to 2^. 

improvement comes from the split marking strategy explicitly taking into account the 
error in the solvent and trying to control where it is large. 

The second plot in Fig. 7 shows that solving the dual problem using piecewise qua- 
dratic elements has substantial additional cost. Although, it is likely that this cost would 
be mitigated by the additional work needed when solving the nonlinear RPBE (see 
Eq. 2.5). In contrast, the goal-oriented strategy using an indicator constructed from a 
linear dual problem (and split marking) does not suffer from the same problem, and is 
the most efficient method when the total cost is taken into account. 

6.2. Performance of Preconditioners. As was discussed in section 5, classical multi- 
grid (MG) preconditioners perform best in the uniform refinement setting, where there is 
a rapid geometric growth in the number of unknowns as the mesh is refined. In the lo- 
cal refinement setting, this growth is much slower, and is frequently subgeometric when 
the refinement is concentrated in the neighborhood of a low dimensional feature (e.g., a 
point or a line). As a result, the per-iteration complexity of classical MG may fail to scale 
linearly (or suffer from a large scaling constant) as the number of unknowns increases. In 
contrast, many local MG preconditioners do not have the same restriction, and maintain 
optimal per-iteration complexity in both the local and uniform refinement settings. 

In Figure 8, the number of primal unknowns is shown as a function of the refinement 
level for three different refinement strategies. For the energy-based marking/refinement 
strategy, the growth in the number of unknowns is geometric, but subuniform. For goal- 
oriented refinement the growth is even slower, but still geometric. The dashed reference 
line shows a geometric growth rate proportional to 2*". 

As a measure of the relative locality of each preconditioner, we compute the ratio of 
the number of unknowns processed by the smoother, shown in Figure 9. Classical MG is 
a global method, so the ratio of smoothed to total unknowns is 1 in both the energy based 
and goal oriented refinement cases. As expected, the BEK preconditioner consistently 
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Figure 9. Ratio of the number smoothed unknowns, Xj, to total un- 
knowns, Nj, on each level for different multilevel preconditioners 




Figure 10. Conjugate gradient method iteration counts for the precon- 
ditioners used. 

smooths more unknowns than HB, but fewer than MG. This is because the set processed 
by the BEK is a superset of the set processed by HB. 

By design, the local multilevel preconditioners considered here maintain an optimal 
per-iteration complexity in both the local and global refinement settings. The primary 
challenge for these preconditioners is achieving a bounded condition number, indepen- 
dent of problem size. It can be shown that for a given tolerance, the number of conjugate 
gradient (CG) iterations can be bounded by a function of the condition number. Hence, 
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if the condition number is bounded, so is the number of CG iterations. In Figure 10, we 
report the number of CG iterations as a function of refinement level for MG, BEK, and 
HB. As predicted by theory, application of the HE preconditioner leads to a slow growth 
in the number of CG iterations as the mesh is refined regardless of indicator type. In con- 
trast, for both goal-oriented and energy-based refinement the iteration count is bounded 
for both the classical MG and BEK preconditioner s. The iteration counts for BEK are 
modestly higher than classical MG, but the work per level is reduced, since BEK smooths 
only a fraction of the unknowns smoothed by classical MG. For this reason, BEK is a 
compelling alternative to classical MG and HB. 

7. Conclusion 

In this article, we developed goal-oriented error indicators for accurate computation 
of the solvation free energy from solutions of the regularized Poisson-Boltzmann equa- 
tion. We found that due to oscillations and imbalanced cancellation in the error contribu- 
tions, global marking strategies based on goal-oriented error indicators were not viable 
for driving adaptive mesh refinement. To address this problem, we developed a split 
marking strategy based on considering each subdomain individually. In numerical ex- 
periments, we calculated the solvation free energy for a 921 -atom Fasciculin-1 protein. 
Through these experiments, we showed that the new marking strategy, combined with 
goal-oriented refinement, is more efficient than energy-based refinement in the context 
of solvation free energy calculations. 

The use of adaptive mesh refinement puts a greater burden on the preconditioner to 
maintain optimal runtime efficiency. To address this issue, we investigated the use of lo- 
cal multigrid methods, which have a lower per-iteration complexity compared to classical 
global multigrid. In particular, the BEK variant proved to be a compelling alternative to 
classical multigrid since it has optimal per-iteration complexity, while still maintaining 
a bounded iteration count as the mesh is refined. The result is an iterative solver with an 
optimal overall complexity, scaling linearly with problem size. 
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